
import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns

import skimage.io as io
from skimage import util

from glob import glob
import trackpy

import ast
from ipywidgets import interact, widgets


# import code from the src
### import code from src folder
import sys
sys.path.append('../src/')
import chow_test_v

import warnings
warnings.filterwarnings('ignore')

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error


from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf


file_list=glob("../data/*_Dark_prepared.tif")
images_dict={}
for file_name in file_list:
    images_prepared=io.imread(file_name)
    file_prefix=os.path.split(file_name)[1].split('.')[0].split("_prepared")[0]
    
    binary_file=file_name.split(".tif")[0]+'_segmented'+'.tif'
    images_binary=io.imread(binary_file)
    tracks=pd.read_csv("../data/tracks_overlay_id_"+file_prefix+'.csv')
    tracks_big=pd.read_csv("../data/tracks_overlay_id_"+file_prefix+'_big.csv')
    img_dict={"img":images_prepared,"binary":images_binary,"tracks":tracks,"tracks_big":tracks_big}
    
    images_dict[file_prefix]=img_dict


images_dict.keys()

dict_keys(['43125_Dark', '43121_Dark', '43122_Dark', '43123_Dark'])


def object_anchor_area(traj,binary_images,ax):
    """ calculate overlapped area i.e sum of pixels that are stay the same across time """
    
    min_row,min_col=np.min(traj.loc[:,'bbox-0':'bbox-1'])
    max_row,max_col=np.max(traj.loc[:,'bbox-2':'bbox-3'])
    
    time_start=np.min(traj.frame)
    time_end=np.max(traj.frame)
    binary_img_region=binary_images[time_start:time_end,min_row:max_row,min_col:max_col]
    
    # calculate anchor mask
    anchor_points=np.zeros_like(binary_img_region[0],dtype=np.int32)
    anchor_points=np.sum(binary_img_region,axis=0)
    anchor_mask=anchor_points>=(255*(binary_img_region.shape[0]))
    ax.imshow(anchor_points,vmin=255,vmax=np.max(anchor_points),cmap='plasma')
    
    ax.axis("off")
    
    # overlap
    sum_overlap_px=np.sum(anchor_mask)
    if np.sum(anchor_mask)>0:
        ax.set_title("file_id={} \n particle_id={} \n anchored=True \n overlapped_area={:.2f}".format(traj["file_name"].iloc[0],traj["particle"].iloc[0],sum_overlap_px/np.mean(traj.area)))
    else:
        ax.set_title("{} \n anchored=False".format(traj["particle"].iloc[0]))
    #print(traj.particle.iloc[0],np.mean(traj.area),sum_overlap_px)
    
    result={'particle':traj.particle.iloc[0],
            'mean_area':np.mean(traj.area),
            'overlapped_area':sum_overlap_px,
            'overlapped_area_ratio':sum_overlap_px/np.mean(traj.area),
            'traj_len':traj.shape[0]}
    
    return result


Tracks_len=5
sns.set(font_scale=2)
Anchor_all=pd.DataFrame()

for file_prefix in images_dict.keys():
    print(file_prefix)
    
    tracks=images_dict[file_prefix]["tracks_big"]
    binary_images=images_dict[file_prefix]["binary"]

    # traj longer then 5 frames
    tracks_length=tracks.groupby('particle').apply(len)
    particle_select=tracks_length[tracks_length>Tracks_len].index.values
    tracks_filters=tracks[tracks['particle'].isin(particle_select)]
    
    Ngroups=tracks.groupby('particle').ngroups
    Nrows=np.int(np.ceil(Ngroups/5))
    fig,ax=plt.subplots(nrows=Nrows,ncols=5,figsize=(40,Nrows*7.5))
    suptitle = plt.suptitle(file_prefix, y=1.02)
    plt.tight_layout()
    axx=ax.ravel()


    counter=0
    Anchor_data=[]
    for id, traj in tracks.groupby('particle'):
        result=object_anchor_area(traj,binary_images,axx[counter])
        counter=counter+1
        Anchor_data.append(result)

    Anchor_df=pd.DataFrame(Anchor_data)
    Anchor_df["type"]="fixed"
    Anchor_df["type"][Anchor_df["overlapped_area_ratio"]<=0.50]="rotating"
    Anchor_df["type"][Anchor_df["overlapped_area_ratio"]==0]="moving"
    Anchor_df['file_prefix']=file_prefix
    Anchor_df.to_csv("../results/Anchor_area_{}.csv".format(file_prefix))
    #plt.savefig("../results/Projection_plot_{}.png".format(file_prefix), bbox_extra_artists=(suptitle,), bbox_inches="tight")
    
    
    Anchor_all=pd.concat([Anchor_df,Anchor_all])

43125_Dark
43121_Dark
43122_Dark
43123_Dark


dxx=1.42 # micron/px

fig,ax=plt.subplots(figsize=(20,8))
Anchor_all["mean_area_micron2"]=Anchor_all["mean_area"]*dxx*dxx
my_order=['43121_Dark', '43122_Dark', '43123_Dark','43125_Dark']
sns.boxplot(x='file_prefix',y="mean_area_micron2",hue='type',data=Anchor_all,order=my_order,showfliers = False,ax=ax)
ax.set_title("Mean Area of the Aggregates")
ax.set_ylabel("mean area in microns^2")
ax.set_xlabel("experiment")
plt.savefig("../results/MeanArea_Aggregate_{}.png",bbox_inches="tight")


Moving_particles=Anchor_all[Anchor_all["type"]=='moving']


Moving_particles=Moving_particles.reset_index()
for file_prefix in my_order:
    selection=Moving_particles[Moving_particles["file_prefix"]==file_prefix]
    tracks=images_dict[file_prefix]["tracks_big"]
    binary_images=images_dict[file_prefix]["binary"]
    
    Nparticles=selection.shape[0]
    Nrows=np.int(np.ceil(Nparticles/5))

    
    fig,ax=plt.subplots(nrows=Nrows,ncols=5,figsize=(40,Nrows*7.5))
    suptitle = plt.suptitle(file_prefix, y=1.02)
    plt.tight_layout()
    
    axx=ax.ravel()
    counter=0
    for k in selection.index:
        particle=selection.loc[k]['particle']
        traj=tracks[(tracks["file_name"]==file_prefix) & (tracks["particle"]==particle)]
        result=object_anchor_area(traj,binary_images,axx[counter])
        counter=counter+1

    plt.savefig("../results/Moving_particles_overlay_{}.png".format(file_prefix),bbox_inches="tight")


tracks[(tracks["file_name"]==file_prefix) & (tracks["particle"]==particle)]


tracks["file_name"].unique()

array(['43123_Dark'], dtype=object)


def plot_trajectory_particle_XY_view(particle,ax):
    traj=tracks[tracks["particle"]==particle]
    traj=traj.sort_values(by=["frame"])
    
    ax.plot(traj['x'],traj['y'],'-o',label="particle={}".format(particle))
    ax.set_title("Trajectory")
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    
    
    # plot anomalies
    anomalies=simple_anomaly_detector(traj,ax,5,1.5,plot_flag=False)
    traj_anomalies=traj[traj.frame.isin(anomalies.index)]
    
    #print("anomalies R","particle={} frames={}".format(particle,anomalies.index),traj_anomalies.x)
    ax.plot(traj_anomalies['x'],traj_anomalies['y'],"ro", markersize=5,label="anomalies")
    ax.legend()
    
def plot_trajectory_particle_R(particle,ax):
    traj=tracks[tracks["particle"]==particle]
    traj=traj.sort_values(by=["frame"])
    start_x=traj['x'].values[0]
    start_y=traj['y'].values[0]
    traj['r']=np.sqrt((traj['x'].values-start_x)**2+(traj['y'].values-start_y)**2)

    ax.plot(traj['frame'],traj['r'],'-o',label="particle={}".format(particle))
    ax.set_title("Radius vector")
    ax.set_xlabel("frame")
    ax.set_ylabel("Radius vector")
    
    # plot anomalies
    anomalies=simple_anomaly_detector(traj,ax,5,1.5,plot_flag=False)
    traj_anomalies=traj[traj.frame.isin(anomalies.index)]
    
    ax.plot(traj_anomalies.frame,traj_anomalies.r,"ro", markersize=5,label="anomalies")
    ax.legend()
    
def plot_trajectory_particle_V(particle,ax):
    traj=tracks[tracks["particle"]==particle]
    traj=traj.sort_values(by=["frame"])
    start_x=traj['x'].values[0]
    start_y=traj['y'].values[0]
    r=np.sqrt((traj['x'].values-start_x)**2+(traj['y'].values-start_y)**2)
    traj['r']=r
    traj["speed"]=np.nan
    traj["speed"][1:]=np.diff(r)
    
#   ax.plot(traj['frame'],traj["speed"],'-o',label="particle={}".format(particle))
    ax.set_title("Speed")
    ax.set_xlabel("frame")
    ax.set_ylabel("Speed")
    
    anomalies=simple_anomaly_detector(traj,ax,5,1.5,plot_flag=True)
    traj_anomalies=traj[traj.frame.isin(anomalies.index)]
    ax.plot(traj_anomalies.frame,traj_anomalies.speed,'ro')
    ax.legend()
    

def simple_anomaly_detector(traj,ax,window=3,scale=1.95,plot_flag=False):
    """
    simple anomaly detector
    """
    break_id=[]
    start_x=traj['x'].values[0]
    start_y=traj['y'].values[0]
    r=np.sqrt((traj['x'].values-start_x)**2+(traj['y'].values-start_y)**2)
    traj['r']=r
    traj["speed"]=np.nan
    traj["speed"][1:]=np.diff(r)
    
    # set minimal size of trajectories for anomaly detection
    anomalies=pd.Series('NA',index=[])
    if (traj.shape[0]>5):
        tt=traj.copy()
        tt=tt.sort_values('frame')
        tt=tt.set_index('frame')
        tt=tt['speed']
    
        
        # find candidates for anomaly
        rolling_mean = tt.rolling(window=5).mean()
        mae=mean_absolute_error(tt[window:], rolling_mean[window:])
        deviation = np.std(tt[window:] - rolling_mean[window:])
        lower_bond = rolling_mean - (mae + scale * deviation)
        upper_bond = rolling_mean + (mae + scale * deviation)
        
        
        anomalies_lower = tt[tt<lower_bond]
        anomalies_upper = tt[tt>upper_bond]
        anomalies=pd.concat([anomalies_upper,anomalies_lower])
        
       
        if (plot_flag==True):
            # plot anomalies
            #fig,ax=plt.subplots(figsize=(8,5))
            ax.set_title("trajectory id = {}".format(traj["particle"].values[0]))
            ax.plot(rolling_mean, "g",lw=3,label="Rolling mean trend") # label="Rolling mean trend"
            ax.plot(tt, color="darkblue", label="data")
            ax.plot(upper_bond, "r--",label="Upper Bond / Lower Bond") #label="Upper Bond / Lower Bond"
            ax.plot(lower_bond, "r--")
#             ax.plot(anomalies_upper, "ro", markersize=5,label="anomalies")
#             ax.plot(anomalies_lower, "ro", markersize=5)
            ax.legend()
    return anomalies


moving_particle=Anchor_df[Anchor_df["anchored"]==0]["particle"].values
sns.set(font_scale=1)

fig,ax=plt.subplots(nrows=moving_particle.shape[0],ncols=3,figsize=(15,moving_particle.shape[0]*5))
suptitle = plt.suptitle(file_prefix, y=0.9)

counter=0
for particle_id in moving_particle:
    plot_trajectory_particle_XY_view(particle_id,ax[counter,0])
    plot_trajectory_particle_R(particle_id,ax[counter,1])
    plot_trajectory_particle_V(particle_id,ax[counter,2])
    #seed_analysis(particle_id,ax)
    counter=counter+1
    
plt.savefig("../results/Moving_trajectories_{}.png".format(file_prefix), bbox_extra_artists=(suptitle,), bbox_inches="tight")


def orientation(traj,ax):
    particle=traj['particle'].values[0]
    time=traj['frame'].values
    
    # elivate big jumps of angle 
    # due to the transition from [pi/2,delta+pi/2] or [-pi/2,-delta-pi/2]
    delta=2
    d_orientation=np.diff(traj.orientation)
    d_orientation=np.where((d_orientation)>delta,d_orientation-np.pi,d_orientation)
    d_orientation=np.where((d_orientation)<-delta,d_orientation+np.pi,d_orientation)
    angle=np.cumsum(d_orientation)
    angle=angle*180/np.pi
    
    #========================#
    # median eccentricity
    eccentricity=np.median(traj['eccentricity'])
    
    # max diff angle
    max_diff_angle=np.max(np.abs(d_orientation))*180/np.pi
    
    # cumsum
    cumsum_angle=angle[-1]
    #========================#
    
    result={"particle":particle,
            "med_eccentricity":eccentricity,
            "max_diff_angle":max_diff_angle,
            "cumsum_angle":cumsum_angle}
    
    #angle=angle.rolling(3).median().values
    s = pd.Series(angle)
    #pd.plotting.autocorrelation_plot(s,ax=ax)
    #ax.legend(str(particle))
    ax.plot(traj.frame[1:],angle,'-o',lw=3,label='particle={:}\ne={:.2f}'.format(particle,np.median(eccentricity)))
    ax.legend()
    
    
    return result


##==============================================================##
sns.set(font_scale=2)
Nrows=np.int(np.ceil(Ngroups/5))
fig,ax=plt.subplots(nrows=Nrows,ncols=5,figsize=(40,Nrows*7.5))
suptitle = plt.suptitle(file_prefix, y=1.02)
plt.tight_layout()
axx=ax.ravel()

counter=0
Orientation=list()
for id, traj in tracks.groupby('particle'):
    result=orientation(traj,axx[counter])
    Orientation.append(result)
    counter=counter+1
    
Orientation=pd.DataFrame(Orientation)

plt.savefig("../results/Angle_{}.png".format(file_prefix), bbox_extra_artists=(suptitle,), bbox_inches="tight")


Orientation[["med_eccentricity",'max_diff_angle','cumsum_angle']].hist(layout=(1,3),figsize=(15,5));


import ipyvolume as ipv
ipv.figure()
#ipv.volshow(X3d.astype(float))
# add scatter plot
x=X3d[:,0]
y=X3d[:,1]
z=X3d[:,2]

ipv.scatter(x = X3d[:,0], y = X3d[:,1], z = X3d[:,2], size=1.5, color = "red")
# ipv.scatter(x = X3d[(cluster==1),0].astype(float), y = X3d[(cluster==1),1].astype(float), z = X3d[(cluster==1),2].astype(float), size=1.5, color = "blue")
ipv.show()

	frame	Unnamed: 0	index	label	y	x	major_axis_length	minor_axis_length	orientation	area	...	max_intensity	mean_intensity	min_intensity	compactness	circularity	contrast	peaks	is_peaks	file_name	particle
9	0	124	124	125	1447.698507	2637.570149	44.003175	19.622468	1.163872	670	...	0.080664	0.052770	0.030342	17.974467	0.699123	0.453330	[[10, 13, 0.07343348]]	True	43121_Dark	11
32	1	304	130	131	1448.539514	2638.141337	42.793624	19.772169	1.083070	658	...	0.080786	0.052659	0.030310	16.992359	0.739531	0.454351	[[9, 9, 0.07828066]]	True	43121_Dark	11
64	2	483	130	131	1447.606432	2637.727412	42.790849	19.631727	1.134981	653	...	0.083023	0.053610	0.030397	16.855226	0.745547	0.463995	[[8, 9, 0.078557245], [11, 19, 0.074834675]]	True	43121_Dark	11
95	3	669	134	135	1444.180030	2636.751891	43.154311	19.772284	1.095776	661	...	0.084501	0.053975	0.030364	17.103173	0.734739	0.471314	[[10, 9, 0.07977986]]	True	43121_Dark	11
128	4	853	136	137	1436.439169	2637.314540	44.320408	19.801536	0.848037	674	...	0.080988	0.052219	0.030344	18.026169	0.697118	0.454885	[[10, 9, 0.08098753], [27, 24, 0.07374989]]	True	43121_Dark	11
160	5	1029	130	131	1433.748191	2636.487699	44.882077	19.814419	1.217332	691	...	0.081372	0.053269	0.030320	17.537372	0.716548	0.457081	[[10, 11, 0.076569945]]	True	43121_Dark	11
190	6	1206	124	125	1429.911278	2627.348872	44.865437	19.105635	-1.377613	665	...	0.083593	0.053932	0.030304	17.678625	0.710823	0.467869	[[11, 11, 0.08359334], [12, 29, 0.07092651]]	True	43121_Dark	11
224	7	1388	137	138	1424.060914	2616.868020	51.218115	22.308433	-1.232191	788	...	0.082013	0.053671	0.030303	21.748566	0.577802	0.460403	[[11, 12, 0.082013026], [8, 23, 0.0779645], [1...	True	43121_Dark	11
257	8	1567	132	133	1427.468354	2607.815752	48.246781	19.501770	-0.997746	711	...	0.080656	0.054439	0.030401	19.350859	0.649396	0.452519	[[18, 12, 0.080655865]]	True	43121_Dark	11
289	9	1745	133	134	1427.504762	2597.572789	52.269516	18.931247	-0.675165	735	...	0.080852	0.052300	0.030338	22.299651	0.563523	0.454305	[[24, 12, 0.08085244]]	True	43121_Dark	11
322	10	1932	137	138	1425.472262	2597.728307	52.716605	17.737922	-0.600109	703	...	0.079584	0.051714	0.030338	24.794365	0.506824	0.448012	[[25, 12, 0.079584084]]	True	43121_Dark	11
355	11	2111	128	129	1422.987688	2597.503420	54.313966	18.694574	-0.533041	731	...	0.077682	0.051586	0.030356	25.529876	0.492222	0.438051	[[27, 11, 0.077681825], [17, 17, 0.07501694]]	True	43121_Dark	11
388	12	2283	126	127	1423.031250	2596.906250	50.393967	18.696645	-0.249153	704	...	0.076549	0.051730	0.030315	20.849959	0.602705	0.432649	[[23, 13, 0.076549046], [35, 7, 0.07624147]]	True	43121_Dark	11
421	13	2458	125	126	1429.298772	2596.619372	49.179304	19.837308	0.124155	733	...	0.080598	0.053534	0.030363	19.402067	0.647682	0.452731	[[24, 10, 0.08059812], [36, 9, 0.07887973]]	True	43121_Dark	11
455	14	2637	126	127	1433.891746	2597.380244	48.166567	20.675286	0.639852	739	...	0.081575	0.051336	0.030302	18.991442	0.661686	0.458302	[[20, 16, 0.08157486], [30, 19, 0.07099368]]	True	43121_Dark	11
488	15	2824	141	142	1436.422572	2595.653543	47.493154	22.963412	1.264990	762	...	0.079663	0.051724	0.030309	19.871711	0.632375	0.448784	[[17, 22, 0.07966266], [19, 31, 0.07864646], [...	True	43121_Dark	11
520	16	2996	127	128	1439.142857	2594.166667	46.561048	22.598789	1.558548	756	...	0.091547	0.054747	0.030346	19.074060	0.658820	0.502092	[[18, 25, 0.091546975]]	True	43121_Dark	11
553	17	3175	131	132	1441.069737	2595.580263	47.847639	22.131698	-1.251605	760	...	0.094852	0.054790	0.030393	20.417529	0.615470	0.514667	[[19, 33, 0.0948522]]	True	43121_Dark	11
582	18	3358	135	136	1442.437419	2594.829677	48.432850	22.320783	-1.147310	775	...	0.097383	0.055602	0.030298	19.944426	0.630069	0.525410	[[20, 33, 0.09738266], [24, 24, 0.09219883]]	True	43121_Dark	11
620	19	3531	127	128	1443.083445	2593.368775	46.695231	22.237376	-0.917124	743	...	0.085416	0.053567	0.030344	20.026639	0.627483	0.475746	[[20, 29, 0.08541554], [27, 21, 0.08304623]]	True	43121_Dark	11

Anchor points¶

Read data¶

Calculate objects anchor area¶

Class of the object does not depend on the¶

Plot patriclular class of the objects¶

Speed¶

Orientation¶